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Abstract 

The use of Laplacian eigenbases has been shown to be fruitful in many 
computer graphics applications. Today, state-of-the-art approaches to 
shape analysis, synthesis, and correspondence rely on these natural har- 
monic bases that allow using classical tools from harmonic analysis on 
manifolds. However, many applications involving multiple shapes are ob- 
stacled by the fact that Laplacian eigenbases computed independently on 
different shapes are often incompatible with each other. In this paper, we 
propose the construction of common approximate eigenbases for multiple 
shapes using approximate joint diagonalization algorithms. We illustrate 
the benefits of the proposed approach on tasks from shape editing, pose 
transfer, correspondence, and similarity. 



1 Introduction 



It is well-established that the eigenfunctions of the Laplace-Beltrami operator 
(manifold harmonics) of a 3D shape modeled as a 2-manifold play the role 
of the Fourier basis in the Euclidean space Tau95 , Lev06 . Methods based 



on the Laplace-Beltrami operator have been used in a wide range of applica- 
tions, including remeshing |Kob97, NISA06 , parametrization |FH05| , compres- 
sion |KG00| , recognition |RWPQ5||RusQ7| , and clustering. Many methods in 
computer graphics and geometry processing draw inspiration from the world 
of physics, finding analogies between physical processes such as heat diffusion 



CL06 or wave propagation and the geometric properties of the shape | SOG09|. 



Several papers have studied consistent discretizations of the Laplace-Beltrami 



operator rPP93 , MDSB03 , WMKG08 



The influential paper of Taubin [Tau95 drew the analogy between the clas 



sical signal processing theory and manifold harmonics, showing that standard 
tools in signal processing such as analysis and synthesis of signals can be carried 



Tip Tip Tip Tip TP? 

■02 % 04 05 06 

»JL n? T? "1 «1 

w w w w w 

Figure 1: Top: Laplace- Beltrami eigenfunctions {4>%} and {^i} of cat and lion shapes: 
besides trivial sign flips in eigenfunctions 2 and 3, higher frequency eigenfunctions 
(4-6) manifest quite different behaviors. Bottom: coupled basis function {(pi}, {ipi} 
obtained using approximate joint diagonalization procedure proposed in this paper 
behave consistently. Hot and cold colors represent positive and negative values, re- 
spectively. 



out on manifolds. This idea was extended in KR05 and later in Lev06,jLZ09 



who showed a practical framework for shape filtering and editing using the man- 
ifold harmonics transform. In OBCS + 12 , the authors proposed a novel repre- 



sentation of correspondences between shapes as linear maps between functional 
spaces on manifolds. In this representation, the Laplace-Beltrami eigenbases of 
the shapes play a crucial role, as they allow to parametrize the linear map as a 
matrix mapping the Fourier coefficients from one shape to another. 

Applications involving multiple shapes rely on the fact that the harmonic 
bases computed on each shape independently are compatible with each other. 
However, this assumption is often unrealistic. First, eigenfunctions are only 
defined up to sign flips for shapes having simple spectra (i.e., no multiplicity 
of eigenvalues). In this case every vector in the eigenspace is called an eigen- 
vector. In the more general case, eigenfunctions corresponding to an eigenvalue 
with non-trivial multiplicity are not defined at all, and one can select an ar- 
bitrary orthonormal basis spanning each such an eigenspace. Second, due to 
numerical instabilities, the ordering of the eigenfunctions, especially those rep- 
resenting higher frequencies, is not repeat able across shapes. Finally, harmonic 
bases computed independently on different shapes can be expected to be rea- 
sonably compatible only when the shapes are approximately isometric, since 
isometries preserve the eigenfunctions of the Laplace-Beltrami operator. When 
this assumption is violated, it is generally impossible to expect that the n-th 
harmonic of one shape will correspond to the n-th harmonic of another shape. 
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These drawbacks limit the use of harmonic bases in simultaneous shape analysis 
and processing to approximately isometric shapes, they do not allow to use high 
frequencies, and usually require some intervention to order the eigenfunctions 
or solve sign ambiguities. 

Contributions. In this paper, we propose a general framework allowing to 
extend the notion of harmonic bases by finding a common (approximate) eigen- 
basis of multiple Laplacians. Numerically, this problem is posed as approximate 
joint diagonalization of several matrices. Such methods have received limited at- 
tention in the numerical mathematics community BGBM93 and have been em- 
ployed mostly in blind source separation applications CS93 , CS96 , Yer02 , Zie05| . 
Philosophically similar ideas have been used in the machine learning community 
for multi-view spectral clustering IKRDI11 . To the best of our knowledge, this 



is the first time they are applied to problems in computer graphics. We show a 
few example of applications of such coupled quasi-harmonic bases in Section 5. 



2 Background 

Let us be given a shape modeled as a compact two-dimensional manifold X. 
Given a smooth scalar field / on X, the negative divergence of the gradient of 
a scalar field, Af = — divgrad/, is called the Laplace- Beltrami operator of / 
and can be considered as a generalization of the standard notion of the Laplace 



operator to manifolds Tau95 LZ09 



Since the Laplace-Beltrami operator is a positive self- adjoint operator, it ad- 
mits an eigendecomposition with non-negative eigenvalues A and corresponding 
orthonormal eigenfunctions </>, 

Act) = M (1) 

where orthonormality is understood in the sense of the local inner product in- 
duced by the Riemannian metric on the manifold. Furthermore, due to the 
assumption that our manifold is compact, the spectrum is discrete, = Ai < 
A2 < • • • . In physics, ([!]) is known as the Helmholtz equation representing the 
spatial component of the wave equation. Thinking of our shape as of a vibrat- 
ing membrane, the eigenfunctions fa can be interpreted as natural vibration 
modes of the membrane, while the A^'s assume the meaning of the correspond- 
ing vibration frequencies. The eigenbasis of the Laplace-Beltrami operator is 
frequently referred to as the harmonic basis of the manifold, and the functions 
fa as manifold harmonics. 

In the discrete setting, we represent the manifold X as a triangular mesh 
built upon the vertex set {xi, . . . ,x n }. A function / on the manifold is rep- 
resented by the vector f = (/(xi), . . . ,/(x n )) T of its samples. A common ap- 
proach to discretizing manifold harmonics is by first constructing a discrete 
Laplace-Beltrami operator on the mesh, represented as an n x n matrix, fol- 
lowed by its eigendecomposition. A popular discretization is the cotangent 
scheme [MDSBQ3] , resulting in the generalized eigenvalue problem W4> = D$A, 
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or equivalent ly, the eigenvalue problem L$ = $A, where L = D 1 W, 



D = diag(si,...,s n )/3, 

_ J (cot(a^) + cot(fiij))/2 i ^ j; 



E 



(2) 
(3) 



Si > denotes the sum of the areas of all triangles sharing the vertex z, and 
C6ij , Pij are the two angles opposite to the edge between vertices z and j in 



the two triangles sharing the edge |WBH + Q7 RWP06 . Here, $ = (0i, . . . 
is the matrix of eigenvectors arranged as columns, discretizing the eigenfunc- 
tions of the Laplace-Beltrami operator at the sampled points. Numerically, the 
aforementioned eigendecomposition problem can be posed as the minimization 



mm 



off ($ T W X $) s.t. $ T D X $ = I, 



(4) 



CS96 . Several 
method in fact 



of the sum of squared off-diagonal elements, off (A) = 
classical algorithms for finding eigenvectors such as the Jacobi 
try to reduce the off-diagonal values in an iterative way. 

Laplace-Beltrami eigenbases are equivalent to Fourier bases on Euclidean 
domains, and allow to represent square-integrable functions on the manifold as 
linear combinations of eigenfunctions, akin to Fourier analysis. In particular, 
solutions of PDEs on non-Euclidean domains can be expressed in the Laplacian 
eigenbasis, giving rise to numerous efficient methods for computing e.g. local 



descriptors based on fundamental solutions of heat and wave equations SOG09 
ASCII , isometric embeddings of shapes BN03 Rus07 , diffusion metrics CLQ6 



shape correspondence and similarity |RWP05||BBK + 10[|OMMG10||DK10 



Unlike the Euclidean plane where the Fourier basis is fixed, in non-Euclidean 
spaces, the harmonic basis depends on the domain on which it is defined. Many 
applications working with several shapes (such as shape matching or pose trans- 
fer) rely on the fact that harmonic basis defined on two or more different shapes 
are consistent and behave in a similar way Lev06 . While experimentally it is 



known that often low-frequency harmonics have similar behavior (finding pro- 
trusions in shapes, a fact often employed for shape segmentation |ReulQ| ), there 
is no theoretical guarantee whatsoever of such behavior. Theoretically, consis- 
tent behavior of eigenfunctions can be guarant eed only in very restrictive set- 
tings of isometric shapes with simple spectrum OBCS+12 . As an illustration, 



we give here two examples of applications in which the assumption of consis- 
tent behavior of basis functions is especially crucial. Additional applications are 

discussed in the Section 5. 

Functional correspondence. Ovsjanikov et al. OBCS + 12| proposed an 



elegant way to avoid direct representation of correspondences as maps between 
shapes using a functional representation. The authors noted that when two 
shapes X and Y are related by a bijective correspondence t : X — » Y, then 
for any real function / : X — >> R, one can construct a corresponding function 
g : Y — ^ R as g = f o£ _1 . In other words, the correspondence t uniquely defines 
a mapping between two function spaces T : T{X, R) —> T{Y, R), where T{X, R) 
denotes the space of real functions on X. Furthermore, such a mapping is linear. 
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Laplacian eigenbases, T(4>i) — E 7 >o c ij^P' 




Common bases, T(<^) = ^j>o c ^'^ 

Figure 2: Matrix C of coefficients expressing a given correspondence between two 
po ses of an ele phant (left) and elephant and horse (right) in functional representation 
of OBCS + 12 . The shapes are shown in the first row, with similar colors representing 
corresponding points. Second and third rows: coefficients defined using Laplacian and 
our coupled bases, respectively. 



Equipping X and Y with harmonic bases, {4>i}i>i and {ipj}j>i, respectively, 
one can represent a function / : X — >• R using the set of (generalized) Fourier 
coefficients {a^}i>i as / = J^>i a i^i- Then, translating the representation into 
the other harmonic basis, one obtains a simple representation of the correspon- 
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dence between the shapes 



(5) 



i,3>l 



where Cij are Fourier coefficients of the basis functions of X expressed in the 
basis of Y", defined as T(^) = J2j>i Cij^j- The correspondence can be thus by 
approximated using k basis functions and encoded by a k x k matrix C = (c^-) 
of these coefficients, referred to as the functional matrix in OBCS + 12 . In this 



representation, the computation of the shape correspondence t : X —> Y is trans- 
lated into a simpler task of finding C from a set of correspondence constraints. 
This matrix has a diagonal structure if the harmonic bases are compatible, an 
assumption crucial for the efficient computation of the correspondence. How- 
ever, the authors report that in practice the elements of C spread off the diagonal 
with the increase of the frequency due to the lack of perfect compatibility of the 
harmonic bases. 

Figure 3: Pose transfer from horse (leftmost) to camel shape (second from left) by 
substituting the first 6 Fourier coefficients in the decomposition of extrinsic coordinates 



of the shape in the Laplacian eigenbasis as done in Lev06 (third from left) and joint 
approximate eigenbasis (rightmost). Coupling between the basis function is crucial for 
this approach to work. 



Pose transfer. Levy Lev06 proposed a pose transfer approach based 
on the Fourier decomposition of the manifold embedding coordinates. Given 
two shapes X and Y embedded in R 3 with the corresponding harmonic bases 
{</>i}i>i and {^i}i>i 7 respectively, one first obtains the Fourier decompositions 
of the embeddings 

X = £>&, Y = £>^ (6) 

i>l i>l 

(we denote by X and Y the Euclidean embeddings of manifolds X and Y, and 
by a^,t>i the three-dimensional vectors of the Fourier coefficients corresponding 
to each embedding coordinate). Next, a new shape Z is composed according to 

n 

j — l i>n 



6 



with the first n low frequency coefficients taken from X, and higher frequencies 
taken from Y. This transfers the "layout" (pose) of the shape X to the shape 
Y while preserving the geometric details of Y. This method works when the 
a^'s and the b^'s are expressed in the same "language", i.e., when the Laplacian 
eigenfunctions behave consistently in X and Y. 



3 Coupled quasi-harmonic bases 

Let us be given two shapes X, Y with the corresponding Laplacians Ax , Ay . Q 
If X and Y are related by an isometry t : X — >• Y and have simple Laplacian 
spectrum (no eigenvalues with multiplicity greater than 1), the eigenfunctions 
are defined up to a sign flip, ipi = ±<pi o If some eigenvalue A* = . . . = Ai +P 
has multiplicity p + 1, the individual eigenvectors are not defined, but rather 
the subspaces they span: span{^, . . . , ^i+p} = span{<^, . . . , 4>i+ p } o t -1 . More 
generally, if the shapes are not isometric, the behavior of their eigenfunctions 
can differ dramatically (Figure [lj top). 

This poses severe limitations on applications we mentioned in the previ- 
ous section. Figure [2] (middle) shows the coefficient matrix C defined in ([5|, 
representing the functional correspondence between two shapes. For near- 
isometric shapes (two deformations of an elephant, Figure [5J middle left), since 
ipi ~ ot _1 , the coefficients Cij ~ ±5^, and thus the matrix C is nearly diag- 
onal. However, when trying to express correspondence between non-isometric 
shapes (elephant and horse, Figure |5J middle right), the Laplacian eigenfunc- 
tions manifest a very different behavior breaking this diagonality. 

The same problem is observed when we try to use the technique of Levy 



Lev06 for pose transfer by expressing the embedding coordinates of the shape 
in the respective Laplacian eigenbasis and substituting the low-frequency coef- 
ficients from another shape. Levy disclaims that his method works "provided 
that the eigenfunctions that correspond to the lower frequencies match" |LevQ6| . 
However, such a consistent behavior is not guaranteed at all; Figure [3] (third 
from left) shows how the pose transfer breaks when the eigenfunction are incon- 
sistent. 

The main idea of this paper is to try to find bases ^ that approximately 
diagonalize the respective Laplacians (Ax4> ~ Ax</>, Ay^ « Ay0) and are 
coupled (ifji « (pi oT 1 ). Such new coupled bases, while being nearly harmonic, 
make the basis functions consistent across shapes and cure the problems we 
outlined above (Figure [l] bottom). Figure [2] (bottom) shows the functional 
correspondence represented in the coupled bases {(f>i}i>i : {i/Ji}i>i- Due to the 
coupling, the basis functions behave consistently resulting in almost perfectly 
diagonal matrices C even when the shapes are highly non-isometric (bottom 
right). Likewise, Figure [3] (rightmost) shows that pose transfer using Fourier 
coefficients in the coupled bases works correctly for shapes with inconsistent 
Laplacian eigenfunctions. 



^■We consider the case of a pair of shapes for the mere sake of simplicity. Extension to a 
collection of more shapes is straightforward. 
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Approximate joint diagonalization. We assume that the shapes are 
sampled at nx,ny points, and their Laplacians are discretized as matrices 
Wj,Wy and Dj , Dy of size tlx x nx and ny x ny respectively, as defined 
m ([3|. We further assume that we know the correspondence between / points 
Xj 1 , . . . , Xj t and yy , . . . , yy (as we show in Section 5, very few and very rough 
correspondences are required). We represent these correspondences by an I x nx 
matrix P with elements pij. = 1 and zero elsewhere; the I x ny matrix Q is 
defined accordingly as q^y = 1 zero elsewhere. 

The problem of joint approximate diagonalization (JD) can be formulated 
as the coupling of two problems Q, 

min off ($ T W X $) + off (# T W Y #) + /i||P$ - Q*||| 

s.t. * T D X *=I, ^ T D Y #=I (8) 

where off denotes some off-diagonality penalty, e.g., the sum of the squared 
off-diagonal elements as defined in Section 2. The parameter \i determines the 
coupling strength. For \i = 0, the problem becomes uncoupled and boils down 
to individual diagonalization Q of the two Laplacians. The joint eigenvectors 
in this setting coincide with the eigenvectors of the Laplacians: $ = $ and 

We can parametrize the joint basis functions as linear combinations of the 
Laplacian eigenvectors, $ = 4>A and \£ = \£B, where A and B are matrices of 
combination coefficients of size nx x nx and ny x ny, respectively. Noticing 

* T ~ T 

that $ Wj$ A T $ T Wj$A A T A X A, and same way, ^ Wyf = B T A y B, 

we transform problem ([8| into 

min off (A T A x A) + off (B T A Y B) + /i||P$A - Q^B||| 

A,B 

s.t. A T A = I, B T B = I (9) 

Since A and B are orthonormal, they act as isometries in the respective eigenspaces 
of the Laplacians of X and Y . We can thus think geometrically of problem Q 
as an attempt to rotate and reflect the eigenbases $ and ^ such that they align 
in the best way (in the least squares sense) at corresponding points, while still 
approximately diagonalizing the Laplacians. 

Since in many applications we are not interested in the entire eigenbasis 
but in the first k eigenvectors, this formulation is especially convenient, as it 
allows us to express the first k joint eigenvectors as a linear combination of k' 
eigenvectors (we provide a justification of this assumption in Appendix A), thus 
having the matrices A and B of size k' x k: 

min off (A T A x A) + off (B T A y B) + /x||P*A - Q*B||| 

A,B 

s.t. A T A = I, B T B = I, (10) 

where $ = {(f>i, . . . ,<j>k') and Ax = diag(A^-, . . . , A^); matrices ^, Ay are de- 
fined accordingly. Typically, k,k' -C nx,ny, and thus the problem is much 
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smaller than the full joint diagonalization (|8J) - Note that the coupling term pro- 
vides Ik constraints, so in order not to over-determine the problem, we should 
have 2k' > I. Typical values used in our experiments were k' ~ 30, I ~ 15, and 
k ~ 20. 



Off-diagonality penalty. It is important to note that in problem (10) 
the use of the sum of squared off-diagonal elements as the off-penalty does not 
produce an ordered set of approximate joint eigenvectors. This issue can be 
solved by using an alternative penalty, 

||* T L X $ - A x ||| = ||A T A X A - A x |||, (11) 

which is similar to the sum of squared off-diagonal elements but also includes 
the difference of the diagonal elements. The penalty for the shape Y is defined 
in the same way. 

Coupling term. It is possible to change the coupling strength at differ- 
ent frequencies, by using a more generic form of the coupling term |(P$A — 
Q\&B)V|||, where V = diag(^i, . . . ,Vk) is a diagonal matrix of weights V{ de- 
creasing with frequency. Such frequency-dependent coupling can be useful in 
applications like pose transfer or functional correspondence we are interested 
in strong coupling at low frequencies and can afford weaker coupling of higher 
ones. 

Procrustes problem. Interestingly, in the limit case \i — >> oo where we 



can ignore the off-diagonal penalties, problem ( 10 ) becomes 



min||P$A-Q#B||! s.t. A A = I, B B = I, (12) 

A,B 

Using the invariance of the Frobenius norm under orthogonal transformation, 



we can rewrite problem ( 12 ) as an orthogonal Procrustes problem 



nun ||P* -Qtf ft ||| s.t. ft T ft = I, (13) 

where ft = BA T . The problem has an analytic solution ft = SR T , where 
* P T Q\& — SER T is the singular value decomposition of the matrix * P T Q\£ 



with left- and right singular vectors S.R Sch66 . Then, A = S and B = R 



4 Numerical computation 

Problem ( [I()| ) is a non-linear optimization problem with orthogonality con- 
straints. In our experiments, we used the first-order constrained minimization 
algorithm implemented in MATLAB Optimization Toolbox. We provide below 
the gradients of our cost function. 

Gradient of the off-diagonality penalty is given by 

V A ^(A t A x A)2 = 4(A X AA T A X A - O o AA X ), (14) 
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where O is a matrix of equal columns o = diag(A T A^A) and o denotes element- 
wise product of matrices (see derivation in Appendix B). Gradient of the alter- 



native penalty (11) is derived similarly as 

V A ||A T A X A-A X ||| = 4(A x AA t A x A- k x AA x ). (15) 

Gradient of the coupling term w.r.t. to A is given by 

V A ||P^A - Q#B||| = 2$ I P T (P#A - Q^B); (16) 

the gradient w.r.t. B is obtained in the same way. 

Initialization. Assuming the Laplacians are nearly jointly diagonalizable 
and have simple spectrum, their joint eigenvectors will be equal to the harmonic 
basis functions up to sign flips. Thus, in this case A and B are diagonal matrices 
of ±1. This was found to be a reasonable initialization to the joint diagonaliza- 
tion procedure. We set A = I, and then solve sign flip by setting the elements 
of B to 

f +1 z = j, ||P^-Q^|| < ||P&+Q^||; 
bij={ -1 i = j, ||P&-Q^|| > ||P&+Q^||; (17) 
[ else. 

Initialized this way, the problem starts with decoupled but diagonalizing bases, 
and the optimization tries to improve the coupling. 

Band-wise computation. In applications requiring the computation of 
many approximate joint eigenvectors (k ^> 1), rather than solving problem ( |To| ) 
for k x k! matrices A,B, we can split the eigenvectors into non-overlapping 
bands of size k" \ and solve k/k" problems ( fToj ) for k" x k" matrices A,B (here 
we assume that A: is a multiple of k" and that k" > I). For the zth band, we use 
$ = (4>(j-i)k", • • -,4>ik") and &x = diag(A^_ 1)A .„, . . . , \f k n) and \&, A Y defined 



accordingly in (10). 



5 Results and Applications 



In this section, we show additional examples of coupled bases construction, as 
well as some potential applications of the proposed approach. We used shapes 
from publicly available datasets |BBKQ8[[SPQ4l|SMKFQ4] . Mesh sizes varied 



widely between 600 - 25K vertices. Discretization of the Laplace-Beltrami op- 
erator was done using the cotangent formula MDSB03 . In all our examples, 



we constructed coupled bases solving the JD problem ( 10 ) with off-diagonality 
penalty ( pTj ), as described in Section 4. Typical time to compute 15 joint eigen- 
vectors was about 1 minute. 

Figure [4] (top) shows examples of joint diagonalization of Laplacians of dif- 
ferent shapes: near isometric (two poses of an elephant) and non-isometric (ele- 
phant and horse). We computed the first k = 20 joint approximate eigenvectors. 
The Laplacians are almost perfectly diagonalized by the obtained coupled bases 
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Figure 4: Examples of joint diagonalization of Laplacians of near- isometric shapes 
(two poses of an elephant, top left) and non-isometric shapes (elephant and horse, 
top right; four humanoids, bottom). Second and fourth rows show the approximate 
diagonalization of the respective Laplacians. Coupling was done using 40 points for 
the elephant and horse shapes and 25 for the humanoid shapes. Numbers show the 
ratio of the norm of the diagonal and off-diagonal values. 



in the case of near-isometric shapes; for non-isometric shapes, off-diagonal ele- 
ments are more prominent. Nevertheless, a clear diagonally-dominant structure 
is present. Figure [4] (bottom) shows an additional example of joint diagonaliza- 
tion of the Laplacian matrices of four humanoid shapes by first 15 elements of 
the coupled bases. 

Sensitivity to correspondence error. Figure [5] shows the sensitivity 
of the presented procedure to errors in correspondence used for coupling. In 
this experiment, we used coupling at 10 points that deviated from groundtruth 
correspondence by up to 15% of the geodesic diameter of the shape. Figure [6] 
shows examples of the first elements of the coupled bases obtained in the lat- 
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Figure 5: Sensitivity of joint diagonalization to errors in correspondence (in % of 
geodesic diameter of the shape) used in the coupling term. Shapes are shown with 
similar colors representing corresponding points. Correspondences between 10 points 
used for coupling are shown with lines. Numbers show the ratio of the norm of the 
diagonal and off- diagonal values. 




Figure 6: Coupled bases elements of the human shapes from Figure [5] obtained us- 
ing 10 points with inaccurate correspondence (error of 15% geodesic diameter) for 
coupling. 



ter setting. This experiment illustrates that very few roughly corresponding 
points are required for the coupling term in our problem, and that the proposed 

procedure is robust to correspondence noise. 

Shape correspondence. Ovsjanikov et al. OBCS + 12| computed cor- 
respondence between near-isometric shapes from a set of constraints on the 
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k x k functional matrix C (encoding the correspondence represented using the 
first k harmonic functions as described in Section 2). Given a set of functions 
/ij • • • 5 fp on X and corresponding functions gi,...,g q on Y (for example, in- 
dicator functions of stable regions) represented as nx x p and ny x p matrices 
F = (fi, . . . ,f p ) and G = (gi, • • • ,gp), respectively, C is recovered from a system 
of pk equations with k 2 variables, 



F T 4> = G T #C. 



(18) 



Additional constrains stemming from the properties of the matrix C are also 
added IOBCS+12 . 



The use of our coupled bases in place of standard Laplace-Beltrami eigen- 
bases allows to exploit the sparse structure of C, which is mentioned by Ovs- 
janikov et al., but never used explicitly in their paper. Approximating C « 
diag(cn, . . . , c/e/e), we can rewrite (18) as a system of pk equations with only k 



variables corresponding to the diagonal of C, 
diag(g^) 



diag(gjtf) 



en 



Ckk 



- T 



- T 



(19) 



Equation (19) allows to use significantly less data to fully determine the cor- 
respondence, and is also more computationally efficient. The method can be 
applied iteratively, by first establishing some (rough) correspondence to pro- 
duce coupled bases, then finding C by solving (19) and using the computed 
correspondence to refine the coupling and recompute the bases. 

Figure [7] shows an example of finding functional correspondence between 
non-isometric shapes of human and gorilla. As functions {/i}f =1 and {^}f =1 
we used bi nary ind icator functions of p = 25 regions detec ted using th e MSER 
algorithm LBBlll. We compared the method described in OBCS + 12 for com- 



puting C by solving the system (18) in the Laplace-Beltrami eigenbases (Fig- 
ure [7| middle) and the diagonal-only approximation ( 19 ) in the coupled bases. In 
both cases, we used k = 20 first basis vectors. Then, point-wise correspondence 



was obtained from C using the ICP-like approach proposed in OBCS + 12 . 

Simultaneous mesh editing. Rong et al. |RCGQ8| proposed an approach 
for mesh editing based on elastic energy minimization. Given a shape with 
embedding coordinates X, the method attempts to find a deformation field d 
producing a new shape X' = X + d, providing a set of user-defined n' anchor 
points for which the displacement is known (w.l.o.g. assuming to be the first n' 
points, di = &[ for i = 1, . . . , ??/), as a solution of the system of equations 



k b Lj 



M T 



M " 




d " 




" 











d' 



(20) 



where M = (I, 0) T is annxn' identity matrix, 7 are unknown Lagrange multipli- 
ers corresponding to the constraints on anchor points, and A^, k c are parameters 



13 




Figure 7: Functional correspondence computation. Left: some of the MSER 
regions. Second and third columns: correspondence and matrix C computed 
using the method of OBCS + 12 in Laplace-Beltrami eigenbases and diagonal- 



only approach in coupled bases, respectively. Fourth column: groundtruth cor- 
respondence matrix C represented in coupled bases. 



trading off between resistance to bending and stretching, respectively RCG08 . 
The system of equations can be expressed in the frequency domain using k <Cn 
first harmonic basis functions, 



(k b L 2 x - k c L x )$ 

M T $ 





a 




" 








d' 



(21) 



where a = <& d are the k Fourier coefficients. The desired deformation field is 
obtained by solving the system of equations for a and tran sforming it to the 
spatial domain d = &a (for details, the reader is referred to |RCGQ8| ). 

Using coupled bases, it is possible to easily extend this approach to simul- 



taneous editing of multiple shapes, solving the system (21) with the coupled 
basis $ in place of and applying the deformation to the second mesh us- 
ing d = Figure [8] exemplifies this idea, showing how a deformation of the 
cat shape is automatically transferred to the lion shape, which accurately and 
naturally repeats the cat pose. 

Shape similarity. The diagonalization quality of the Laplacians can be 
used as a criterion for shape similarity, with isometric shapes resulting in ideal 
diagonalization. With this approach, it is possible to compare two shapes from 
a small number of inaccurate correspondences provided for coupling in the joint 
diagonalization problem. It is possible to compare more than two shapes at once, 
with the diagonalization quality acting like some "variance" of the collection of 
shapes. 

Figure [9] shows the similarity matrix between 25 shapes belonging to 8 dif- 
ferent classes. Each shape is present with 3-4 near-isometric deformation. We 
used 25 points for coupling; dissimilarity of a pair of shapes was computed by 
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Figure 8: Simultaneous shape editing in the frequency domain using the approach 
of RCG08 . Top: editing the cat shape (green points denote the anchor vertices used 



as constraints in problem (20)). Bottom: the same pose is automatically transferred 
to the lion shape using coupled basis. 



jointly diagonalizing the respective Laplacians and then computing the average 
ratio of the norms of the diagonal and off-diagonal elements of both matrices. 



6 Conclusions 

We showed several formulations of numerically efficient joint approximate diag- 
onalization algorithms for the construction of coupled bases of the Laplacians of 
multiple shapes. Such quasi-harmonic bases allow to extend many shape anal- 
ysis and synthesis tasks to cases where the standard harmonic bases computed 
on each shape separately cease being compatible. The proposed construction 
can be used as an alternative to the standard harmonic bases. A particularly 
promising direction is non-rigid shape matching in the functional correspon- 
dence representation. The sparse structure of the functional matrix C has not 
been yet taken advantage of for the computation of correspondence in a proper 
way, and naturally calls for sparse modeling methods successfully used in signal 
processing. However, the correctness of this model largely depends on the basis 
in which C is represented, and standard Laplace-Beltrami eigenbasis performs 
quite poorly when one deviates from the isometry assumptions. In follow-up 
studies, we intend to explore the relation between sparse models and joint ap- 
proximate diagonalization in shape correspondence problems. 



Appendix A - Perturbation analysis of JD 

We analyze here the solution of the basic problem (J8|. For simplicity of anal- 
ysis, we assume that nx = rty = Z, the points are ordered as ik = jk — and 
/i — )• oo. This setting of the joint diagonalization problem has been considered 
in numerical mathematics BGBM93 and blind source separation |CS96| liter- 
ature. In this case, the matrices are of equal size, have row- and column- wise 
correspondence, and a single basis (<£ = \£) is searched. 
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Figure 9: Shape similarity using joint diagonalization. Darker colors represent more 
similar shapes. One can clearly distinguish blocks of isometric shapes. Also, two 
classes of two- and four-legged shapes (marked with green and blue) are visible. Small 
figures show representative shapes from each class. 



We further assume that Lx = 4>A<I> T has a simple r-separated spectrum 
(i.e., \Xi — Xj \ > r). If Y is a near-isometric deformation of X, its Laplacian can 
be described as a perturbation Ly = $A4> T + eR of Lx . Ignoring permutation 
of eigenfunctions and sign flips, the joint approximate eigenbasis can be written 
as the first-order perturbation |Car95| 



4>i 



<§>i 



(22) 



where = <pjK<pj /2(Aj 
norm \a id \ < ||R|| 2 /2r = 



— Xi). We can bound these coefficients by the spectral 
Q! max . We can conclude that the first k joint approxi- 
mate eigenvectors <$> x , . . . , (j> k can be well represented as linear combinations of 
01, . . . with square error bounded by c^ j>k , \aij\ 2 < e(n x - k f - l)a max . 
This result also justifies the use of band-wise computation discussed in Section 
4. 

To relate this bound to the geometry of the shape, let us assume that X 
and Y have the same connectivity and that the angles /?, /3 of triangles in the 
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two meshes satisfy Oo < /?, f) < tt — Oo (all triangles are at least Oo fat), the area 
elements are at least s, s > so, and each vertex is connected to at most v vertices. 
We assume that the mesh Y is obtained as a deformation of mesh Y changing 
the angles by \/3 — /3\ < AO and the area elements by 1 — S < s/s < 1 + 5. We are 
interested in a bound on the spectral norm ||Lx — Ly H2 = 6 I|R||2 expressed in 
terms of parameters AO, S (strength of non-isometric deformation) and constants 
#0, so- From trigonometric identities, we have for i ^ j 



Wij\ < 



sin(/3 - f3) 



sin f3 sin f3 



sin AO AO 

< — 2 — < — o — • 
sin Oq sin Oq 



(23) 



and \wu - w u \ < 

\hj ~ hj I = ^IwijSi 
in the bounds on \vbj 



sin 2 0q Applying the triangle inequality, for i / j we get 



_1 | < 3\wij - w ij \s i 1 + 3\wi 
and Si, we get 



< 



< 



w ij s i 
~ W ij\ 

-( — 

so \sin 2 u 
3(A0 + 5) 



Plugging 



AO 



+ cot(# )<H < 



so 



AO 
sin 2 Oq 



+ 



so sin 2 Oo 

from which it follows by norm inequality 



||Lx-L Y || 2 <n 



Lx-L 



Y 1 



< 



6vn 



3/2 



-(AO + 6), 



where e = AO 



So sin c/q 

S is the degree of shape deformation ("lack of isometry"). 



Appendix B - Off-diagonality penalty gradient 

Let us rewrite the penalty as off (A) = ||A T AxA||p - ^(A T AxA)|. The first 
term is bi-quadratic and its gradient derivation is trivial, so here we derive the 
gradient of the second term only. We have J^ i (A T AxA) 2 i = (^2 k a1i X k) ; 
differentiating in a coordinate- wise manner, 



Y ( Y a ki Xk J = 2 Y ( Y a ki Xk J 2a P<l X P S < 



da pq 



iq 



k 

Observing that AA = (a pq X p ) and denoting by o pq = (^Z k al q Xk) the elements 
of the equal-columns matrix O, we get the expression 40 o AA. 
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